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The purpose of this work is to use level set topology optimization to improve the design of 
a representative wing box structure for the NASA common research model. The objective is 
to minimize the total compliance of the structure under aerodynamic and body force 
loading, where the aerodynamic loading is coupled to the structural deformation. A taxi 
bump case was also considered, where only body force loads were applied. The trim 
condition that aerodynamic lift must balance the total weight of the aircraft is enforced by 
allowing the root angle of attack to change. The level set optimization method is 
implemented on an unstructured three-dimensional grid, so that the method can optimize a 
wing box with arbitrary geometry. Fast matching and upwind schemes are developed for an 
unstructured grid, which make the level set method robust and efficient. The adjoint method 
is used to obtain the coupled shape sensitivities required to perform aerostructural 
optimization of the wing box structure. 
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Nomenclature 

sensitivity factor for angle of attack 

vector of Doublet Lattice Method (DLM) box areas 

vector defining influence of wing deformed shape on lift 

compliance of the structure 

pressure coefficient vector 

aerodynamic influence coefficient matrix 

material property tensor 

number of elements attached to a node 

aerodynamic load vector 

body force load vector 

total load vector 

acceleration due to gravity 

element edge length 

indices 

global structural stiffness matrix 

stiffness matrix of an element cut by the boundary 

stiffness matrix of a finite element 

iteration number 
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total lift force 

lift force from built-in twist and camber 

lift force from unit angle of attack 

load factor 

unit normal vector 

adjoint state vector 

dynamic pressure 

aerodynamic stiffness matrix 

force transfer matrix 

displacement transfer matrix 

fictitious time variable 

displacement field or vector 

velocity function 

virtual displacement 

wing box weight 

fixed aircraft weight 

downwash dependent on deformed wing shape 
constant downwash from built-in camber 
point in the design domain 
column vector of V s 
angle of attack 

volume of a cut element that lies inside the structure 
volume of an element 
small number 
structural boundary 

part of boundary subject to displacement boundary conditions 
part of boundary subject to aerodynamic loads 

part of boundary free from boundary conditions and aerodynamic loads 
time step 
strain tensor 

arbitrary vector (shape derivative auxiliary variable) 

material density 

implicit function 

design domain 

structural domain 


I. Introduction 

T opology optimization is the most general form of structural optimization, as it is the least restricted by pre- 
determined design variables. 1 Topology optimization has been used as a tool for improving aircraft wing 
designs. However, the full potential of topology optimization is often not realized, as the design space is limited to a 
pre-determined layout or applied only to local regions or sub-components. For example, Balabanov and Haftka 2 
used truss topology optimization to design the internal structure of a wing, where the discrete ground structure limits 
the design space. The Simple Isotropic Material with Penalization (SIMP) method has been used to optimize a rib 
and spar type wing structure, where the layout of ribs and spars was fixed. 3 The bubble method has been used to 
optimize the layout of a single rib for a set of loadings obtained from pull-up maneuvers and tank pressures. 4 The 
SIMP method has also been used to optimize individual ribs. 5 The spectral level set method has been applied to 
optimize the reinforcement of a wing upper skin to maximize aileron reversal dynamic pressure. 6 Topology 
optimization has also been applied to the design of plate-like wings, through varying the thickness of the plate 7 ’ 8 or 
varying the distribution of different materials. 9 The motivation for this work is to consider the entire 3D wing box 
space as the design domain for topology optimization to explore novel configurations for the entire wing box 
structure. 

Another important consideration in topology optimization of aircraft wings is the aero structural coupling 
between the aerodynamics and the deformed shape of the wing. Early approaches computed the pressure loading for 
an assumed wing shape and did not account for the coupling. 2,4 However, this coupling effect can strongly influence 
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the optimal design. 3 * James 10 used a multidisciplinary framework that optimized both the aerodynamic shape of the 
wing and the topology of the wing box structure using the SIMP method. However, the optimal designs contained a 
significant amount of intermediate density structure, which can be difficult to interpret as a manufacturable design, 
which was also noticed in Ref. 3. An alternative is to use the level set topology optimization method, 11 ' 13 which 
defines the boundary of the structure using an implicit function and produces solutions with clear boundaries. 

Various approaches to level set based topology optimization have emerged in the last decade and there is no one 
method that can be called the level set method in the context of topology optimization. This review intends to 
provide an overview of level set based optimization. More comprehensive reviews can be found in Refs 14-16. 
Common to all level set methods is the representation of the boundary as the zero level set of an implicit function. 
Sensitivity information is then used to update the implicit function and hence the position of the boundary. The 
implicit boundary representation allows natural handling of topological changes, such as merging or insertion of 
holes, although special techniques are usually required to automatically create new holes during optimization. 11 ' 16 
The initial approaches to level set based optimization used shape sensitivity analysis coupled with the level set 
method developed by Osher and Sethian. 17,18 In this approach the implicit function is a signed distance function and 
it is updated by solving a Hamilton- Jacobi type equation using an explicit upwind scheme. 11,12 This is sometimes 
called the conventional or direct approach. 14 Variations on the conventional approach include solving the Hamilton- 
Jacobi equation using semi-implicit additive operator splitting 19 or a semi-Lagrange approach. 20 These methods can 
overcome the Courant-Friedrichs-Lewy (CFL) stability condition that limits the efficiency of the conventional 
approach. Further deviations from the conventional approach use representations of the implicit function other than 
the signed distance function. These include using radial basis functions, which convert the Hamilton- Jacobi PDE 
into an ODE 21 or a simple set of algebraic equations, 22 the spectral level set method, which represents the implicit 
function as coefficients of a Fourier series and methods based on phase fields, ’ where the implicit function is 
used to define the material phase of a point and the update is performed by solving a reaction-diffusion equation. 

The application of the level set method for wing structure optimization has been limited to local regions 6 or did 
not account for the aero structural coupling. 26 In our previous work, the level set method was used to optimize the 
internal structure of a 3D wing box and included the effect of aero structural interaction through coupled analysis and 
coupled sensitivities. 27 However, the effect of body force loading was not included, the level set method was 
implemented on a uniform Cartesian grid and the examples were limited to cuboid design domains. Also, the angle 
of attack was fixed and the trim condition that lift must equal to weight was imposed as a constraint. 

To optimize the internal structure of more realistic wings, such as the Common Research Model (CRM), 28 a 
level set method for more general design domains that are non-rectangular needs to be developed. One approach is 
to immerse the general design domain in a larger uniform level set grid and maintain the shape of the actual design 
domain by fixing the implicit function in the region of its boundary. 29 Another approach is to map the general design 
domain to a uniform grid, for example by using an isoparametric transform. 26 However, this is only applicable when 
a mapping is available and the grid is structured. A third approach is to implement a level set framework that can 
naturally handle unstructured grids within arbitrary domain shapes, which, to our knowledge, has not been done for 
level set based topology optimization. In this work methods are developed to adopt the third approach. 

This paper presents a strategy for optimizing the internal structure of a wing box using the conventional level 
set topology optimization method. The effect of the coupled aerodynamic loading is taken into account through a 
multidisciplinary analysis of the system and a coupled adjoint sensitivity analysis. The design space for the 3D level 
set optimization is the volume contained within the wing box and the optimization can choose any structure within 
that space. Thus, the approach is not limited by any pre-determined structural layout, such as ribs and spars. The 
design space is discretized using a grid fitted to the design domain, which will be non-uniform and may also be 
unstructured. Therefore, level set strategies for non-uniform and unstructured grids are employed. 

II. Level Set Topology Optimization 

A 3D level set method was used to optimize the internal structure of the wing box. A conventional level set 
method was employed. The basic method is briefly outlined followed by details of the implementation for an 
unstructured grid. First, the boundary of the structure is defined as the zero level set of an implicit function: 

0(v) > 0,v E Q 

<0(x) = O,xEr (1) 

0(v) < 0,v ^ Q 
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where (p(x ) is the implicit function and v E Q d , where Q d is the design domain containing the structure, Q C Q d . The 
implicit description of the structure allows the boundary to break and merge naturally, allowing topological changes 
to occur during optimization. 

The level set description of the structural boundary can be used for optimization. 11 ' 13 First, a velocity field is 
derived from shape sensitivity analysis. The implicit function is then updated by solving a Hamilton- Jacobi type 
equation: 


^^ + V0(x,r) — = 0 (2) 

dt v ’ dt 

where t is a fictitious time domain. Equation (2) can be discretized and rearranged to provide a convenient update 
formula for optimization: 


0f +1 =^-Ar|v0f|y nji (3) 

where i is a discrete point within the domain, k is the current iteration number, V n is a velocity function defined 
normal to the boundary, such that a positive velocity moves the boundary inwards, and At is a discrete time step. 
Note that Eq. (3) is the same as solving Eq. (2) using an explicit forward Euler scheme. 18 

A. Fixed Grid Finite Element Analysis 

The design domain is discretized using isoparametric hexahedral elements with the implicit function defined at 
the nodes. The same grid is also utilized as the mesh for the finite element analysis (FEA) of the wing box. Elements 
with incompatible modes are used to analyze the wing box structure, as these elements perform better in bending 
compared with standard tri-linear isoparametric elements. 30 

The implicit function is interpolated within elements using standard tri-linear shape functions, which is useful for 
finding the exact position of the boundary. The implicit function is initialized as a signed distance function, where 
the value of the implicit function indicates the distance to the boundary and the sign is defined by Eq. (1). Details of 
how the signed distance function is computed on an unstructured grid are included in Section II.B. 

For efficiency, the finite element mesh remains fixed during optimization. Therefore, some elements can be cut 
by the structure boundary leading to a discontinuity in material properties within the element. Discontinuous 
elements are treated using a volume-fraction weighted approximation for the element stiffness matrix: 


K c - (fic I P e)K e (4) 

where K c is the approximated stiffness matrix of the cut element, K E is the stiffness matrix of the same element 
without a cut, whose volume is /3 E , and /3 C is the volume of the cut element that lies inside the structure. To compute 
the internal volume of a cut element, /3 C , the hexahedral element is 
first divided into 24 sub-tetrahedra, Fig. 1. The implicit function is 
then interpolated onto the vertices of each sub-tetrahedron. The 
internal volume for each sub-tetrahedron can then be analytically 
computed and summed to obtain total value for the element. This 
is similar to the method proposed by Min and Gibou, 31 although 
their method used five sub-tetrahedra. The choice of 24 sub- 
tetrahedra eliminates directional bias from the interpolation and 
leads to a more robust method, at additional computational 
expense. The method of subdividing elements into simplexes 
(triangles in 2D and tetrahedra in 3D) can also be used to perform 
surface integration when the boundary is represented by an 
implicit function. 30 

The total volume of the element, /5 E , can be efficiently computed using an analytical formula. 31 Elements that lie 
completely outside the boundary are assigned a small volume fraction, 10' 6 . This is to avoid numerical difficulties 
involved when part of the structure becomes completely detached from the main structure, such that it can display 
rigid body motions. 



Figure 1. Dividing a hexahedral element 
into 24 sub-tetrahedra. 
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It is known that using the volume fraction method of Eq. (4) (also called the ersatz material approach) leads to 
reduced accuracy at the boundary due to the approximation of the element stiffness. 33 ' 35 This is a concern for 
conventional level set methods, as the optimization is driven by shape sensitivities defined at the boundary. To 
overcome this problem, sensitivity values are computed at element Gauss points and interpolated at the boundary 
point using a weighted least squares method. 34 The original method uses all Gauss point sensitivities within a 
specified radius of a boundary point. This is modified for a non-Cartesian mesh by using Gauss point sensitivities 
inside elements that neighbor the element where the boundary point is located. A 2D example is shown in Fig. 2. 



Figure 2. Boundary point sensitivity calculation for 2D unstructured mesh. 


Boundary shape sensitivity values are computed at points where the boundary intersects an element edge and at 
grid nodes where the implicit function is zero. Boundary point shape sensitivities are then used to compute the 
velocity function. The derivation of the velocity function from shape sensitivities for the problem studied in this 
work is detailed in Section III.C. 

B. Velocity Extension and Reinitialization 

A key aspect of conventional level set based optimization is the definition of the velocity function from the 
shape derivative. However, the resulting velocity function is only defined on the structural boundary. To update the 
level set function using Eq. (3) discrete velocity values, V n>i are required at all grid points i. For stability, it is 
desirable to define these velocity values such that the signed distance property of the implicit function is maintained 
during the update. 17 ' 18 This can be achieved by extending the velocity to be constant along a line normal to the 
boundary. The velocity extension can be obtained by numerically solving: 18 

(dV n / dt ) + Signup <p/\ |V0|) • W„ = 0 (5) 

One method to solve Eq. (5) is to choose an initial solution V n and then employ an iterative procedure until 
convergence. This is performed in a pseudo-time domain that is different from the time axis used in the main level 
set evolution, Eq. (2). However, this approach can be computationally expensive. On Cartesian or rectilinear grids 
the velocity extension can be more efficiently computed using the fast marching method. 17 The fast marching 
method has also been generalized by Sethian and Vladimirsk 36 to solve Eikonal equations on unstructured grids. 
Eikonal equations take the form: 


Jv0|| = F(jc) in Q (6) 

where F(x) is an arbitrary function. If F(x) = 1, then the solution of the Eikonal equation can be used to reinitialize 
the implicit function as a signed distance function, where the magnitude of the implicit function is found from 
solving Eq. (6) and the sign is defined by Eq. (1). However, the unstructured grid fast marching method for solving 
Eikonal equations requires modification to obtain extension velocities. The method used in Ref. 36 to solve the 
Eikonal equation for reinitialization is first summarized followed by details of how it is modified in this work to 
compute the extension velocities. 
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The fast marching method solves the Eikonal equation using an upwind finite difference scheme. The term 
upwind means that implicit function values are computed from information at points closer to the boundary. Thus, 
the solution is constructed in a progressive manner by marching out from the boundary. Finite differences are used 
to compute the spatial gradients of the implicit function. During reinitialization, the implicit function value at a point 
is computed by solving: 


V0,- • V <j>i = 1 


(7) 


where the information used to construct the gradients is obtained 
from local points closer to the boundary. Thus, these points should 
already have implicit function values computed earlier, as they are 
closer to, or on the boundary. In 3D, up to three local points are 
used, as there are three spatial dimensions. If a grid node is near the 
boundary then the local points are locations where the boundary 
crosses an element edge attached to the node. Interpolation 
between two neighboring grid nodes with opposite implicit 
function sign is used to compute the location of a boundary point. 

If a grid node is not near the boundary, then the local points will be 
the grid nodes connected to the current node that are closer to the 
boundary. Due to the upwind scheme, the local points will have 
known implicit function values, computed during a previous step of 
the fast marching method. 

To solve Eq. (7), a simplex is formed using the current grid 
point and the local points closer to the boundary, Fig. 3. The spatial gradients within the simplex are assumed 
constant and computed using linear shape functions, which is the same as using first order finite differences. For 
example, if the simplex is a tetrahedron, Fig. 3, the gradients of the implicit function are computed by obtaining the 
isoparametric gradient matrix, [B], using the four-node tetrahedron shape functions: 30 



Figure 3. Local calculation of extension 
velocity in 3D. 


V0/-[B] M 


( 8 ) 


where {(p) is a vector of nodal values associated with the tetrahedron, Fig. 3. Using Eq. (8), Eq. (7) becomes: 

V0,.-V0,={0} r [Bf[BM = 1 (9) 

The only unknown in Eq. (9) is the implicit function value at the current grid point, which can be obtained by 
solving a quadratic equation. 

Equipped with a procedure for computing local solutions to Eq. (7), the fast marching method can be 
summarized as: 36 

1. Compute local solutions for grid nodes with local points on the boundary, tag these as “known.” 

2. Find all “unknown” grid nodes connected to the “known” set, tag these as “considered.” 

3. Compute local solutions for the “considered” set of nodes. 

4. Find the node in “considered” with the smallest magnitude of (p , tag this as “trial” 

5. Find all nodes connected to “trial” that are not in “known” and update their local solutions, add these to 
the “considered” set and move the “trial” node to the “known” set. 

6. If the “considered” set is empty, stop. Otherwise return to step 4. 

The algorithm employed to compute extension velocities is developed from the fast marching method used to 
solve the Eikonal equation for reinitialization. The algorithm is similar and employs the same upwind and finite 
differencing techniques. The main difference is that instead of locally solving Eq. (7) for the implicit function value, 
the following is solved for the velocity function value at the current grid node: 

V0,-W„,,=O (10) 
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where the spatial gradient of the velocity function is computed using linear shape functions as before: 

v ^, = [£]R} (id 

where { V n ] is a vector of nodal velocity values associated with the locally constructed simplex. From Eqs. (8 ,10 & 
11) we obtain: 


V0, • W„, = {t\ T [B] T [B]{V n } = 0 (12) 

Due to the upwind scheme, the only unknown in Eq. (12) is the velocity value at the current node. Thus, Eq.(12) can 
be rearranged and solved for the velocity value at the current node. 

Care must be taken when using the fast marching method on an unstructured grid, as the causality in the upwind 
scheme can be invalidated due to the presence of obtuse angles, which can cause numerical instability. 36 In our 
implementation, this problem is partially alleviated by removing local points from the gradient approximation that 
cause large angles at the current point. This is a simple method and is not considered robust, as the resulting spatial 
gradient can be very different from the true gradient. 36 However, for the grids used in this work, experience showed 
this method to be sufficient. More robust procedures can be found in Ref. 36. 

The velocity extension method aims to preserve the signed distance function. However, the implicit function 
may still lose the signed distance property during optimization, due to the approximations in the numerical 
implementation, such as the use of first order finite differences and the heuristic method to deal with obtuse angles. 
Therefore, the implicit function is reinitialized as a signed distance function occasionally during the optimization. 

C. Gradient of the Implicit Function 

The spatial gradient of the implicit function at each grid node is also required to update the structure using Eq. 
(3). Upwind schemes are favored when solving the Hamilton- Jacobi equation, as gradients are computed according 
to the direction that information is traveling, which promotes stability in the solution. 17 ' 18 On a rectilinear grid, 
upwind schemes require a simple choice between forward or backward finite differences. However, on an 
unstructured grid, this is not always possible. 

For an unstructured grid, spatial gradients can be computed using element shape functions in the same way that 
strain values are computed from nodal displacement values in FEA. 30 However, gradients computed using trilinear 
shape functions are only guaranteed to be continuous within an element and can be discontinuous across element 
boundaries. Therefore, as grid nodes are attached to several elements, the gradient calculated at the same node using 
different elements is not certain to be the same. In FEA, a simple remedy is to average the gradient values from each 
element at the node. However, we would like to compute gradients according to an upwind scheme to promote 
stability in the level set method. Therefore, nodal gradients are computed using only the attached element from 
where information flows into the node. To decide which element to use, first the implicit function value at each 
element center is computed. Then, for node i, the following criterion is used to decide which element is used to 
compute the gradient at the node: 


max(sign(V ni )<pj), j = l -e (13) 

where j is an element attached to node i, is the implicit function value at the center of element j and e is the 
number of attached elements. The spatial gradient of the implicit function is then computed at the node using 
information from the chosen element, i.e. the one that satisfies Eq. (13). 

D. The Level Set Topology Optimization Algorithm 

The time step in Eq. (3) is computed using the well-known CFL condition for stability: 

At = min(0.9/j m | n ,/V n; ) (14) 

l 


where h min>i is the shortest element edge length connected to node i. 

A simple heuristic termination criterion is used, where the optimization is terminated if the maximum relative 
change in the objective over the previous ten iterations is less than a small number: 
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(15) 


(c ax - CL ) (cL + CL )<y, ie[k- 9, k] 

where C is the total compliance value, or objective, k is the iteration number and y is a small number. 

Finally, the optimization algorithm is composed of the following major steps: 

1. Assemble the global stiffness matrix from volume ratio weighted element matrices, Eq. (4). 

2. Solve the primary and adjoint equations (Section III.D). 

3. Check for convergence using Eq. (15). Otherwise continue. 

4. Compute shape sensitivities at the boundary (Section II. A) 

5. Velocity extension using the fast marching method (Section II.B) 

6. Spatial gradient calculation using an upwind scheme (Section II. C) 

7. Update the implicit function, Eq. (3) 

III. Aerostructural Optimization Problem 

The objective studied in this work is to minimize the compliance of the internal wing box structure. The loading 
on the structure is composed of the lift force generated by a single flight condition and the weight of the wing. The 
lift distribution is dependent on the deformed shape of the wing, which is in turn dependent on the applied loading. 
Therefore, the structure is part of a coupled, multidisciplinary system. The wing is trimmed by allowing the root 
angle of attack to change so that the total lift force is equal to the total weight of the aircraft. The problem can be 
stated as: 


Minimize : C(u) = fj u 

Subject to : Ku = f t = f a ( u) + N • f g (16) 

L(u) = N-(W b +W c ) 

where u is the deformation vector of the wing, C(u) is the total compliance of the wing, K is the structural stiffness 
matrix, f a {u) is a load vector resulting from the aerodynamic lift force, f g is the body force loading due to the weight 
of the wing, A is a load factor (N = 1 for steady level flight), L(u) is the total lift generated by the deformed wing 
shape, W b is the weight of the wing box that can change during the optimization (and is dependent on the position of 
the structural boundary), whereas W c is the remaining “fixed” weight of the aircraft. The first constraint in Eq. (16) 
is the static equilibrium of the aerostructural system and the second constraint is the trim condition. The 
aerodynamics and pitch-based trim conditions from a tail are not considered for the wing optimization. For the taxi- 
bump load case, the aerodynamic forces are zero and the trim condition is not enforced. 

A. Aerodynamic Modeling 

The subsonic compressible aerodynamics are modeled using the Doublet Lattice Method (DLM). 37-38 The wing 
planform is divided into discrete boxes, each with a control point at the % chord position. The DLM is used to 
compute the pressure acting on each box for a given distribution of downwash at the control points: 


Dc p =w c -a-z + w(u) (17) 

where D is a matrix of non-dimensional aerodynamic influence coefficients and c p is a vector of differential pressure 
coefficients for each DLM box, between the upper and lower wing surfaces. There are three sources of downwash at 
each control point, w c is the slope of the aerofoil centerline (from built-in camber and twist), a is the angle of attack 
at the wing root, z is a column vector of 1 ’ s and w(u) is the additional downwash dependent on the deformed shape 
of the wing. The total lift force generated by the wing is computed from the DLM box planform areas, the pressure 
coefficients and the dynamic pressure: 

L(u) = q • a T c p = q • a T D~ l {w c -a m z + w(u )} (18) 

where q is the dynamic pressure and a is a column vector containing the DLM box areas. 
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B. Coupling the Structural and Aerodynamic Analyses 

The Finite Plate Spline (FPS) method 39 is used to communicate information between the FEA and the DLM 
meshes. A flat plate finite element mesh is constructed between the DLM and FEA discretizations and displacement 
constraints are enforced at locations that coincide with the DLM box centers and a subset of the FEA nodes. Solving 
the FPS finite element equation under these constraints allows the formulation of two transfer matrices that relate 
aerodynamic and structural quantities: 


II 

& 

7s 

(19) 

w(u) = Tu 

(20) 


where S is a matrix that computes the work-equivalent loading on the structure from the pressure coefficients and T 
is a matrix that converts the structural displacement vector into the downwash component dependent on the 
deformed wing shape, also in a work-equivalent way. Discrete Kirchoff Triangular (DKT) plate elements 30 are used 
to form the FPS mesh. 

Combining the trim condition from Eq. (16) with Eqs. (18) and (20) results in an expression for the required 
angle of attack for trim in terms of the deformed shape of the wing and current wing box weight: 

a(u) = [L c +b T u-N-(W b +W c ))/L a (21) 


where L c and L a are constants defined as: 


L c = q • o T D 1 w c 

(22) 

L a = q • a T D~ l z 

(23) 


and b T is a constant vector defined as: 

b T = q-a T D~ l T (24) 

The matrices S and T allow information to pass between the DLM and FEA and they can also be used to 
formulate an aerodynamic stiffness matrix: 


Q = qSD~ l T (25) 

Note that the Q matrix is constant for a given wing planform and flight condition, defined by velocity, air density 
and Mach number. However, the Mach number dependence is contained within D, whereas velocity and air density 
are used to define q. Thus, if the Mach number is changed, then the D matrix has to be recomputed and inverted to 
update the Q matrix. However, if dynamic pressure is changed then the update of Q only involves a scalar 
multiplication of Q. The aerodynamic load vector in Eq. (16) can be computed from Eqs. (17), (19) and (25): 

fa («) = fc ~ a (u) ' f a + Qu (26) 


where a(u) is defined by Eq. (21), f c and /„ are constant load vectors defined as: 


\ 

II 

(27) 

f a =q-SD~ l z 

(28) 


The trim condition in Eq. (16) can now be eliminated and the problem re- stated as: 
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Minimize : C[u) = fj u 
Subject to: Ku = f t = f a (u) + N • f g 


( 29 ) 


where the aerodynamic load vector f a (u) is defined by Eq. (26) and the required angle of attack for trim is defined by 
Eq. (21). The only constraint in Eq. (29) is the static equilibrium equation for the coupled aero structural system, 
which is implicitly satisfied when solving for u. Thus, the optimization problem defined in Eq. (29) can be solved 
using unconstrained methods, although it may be necessary to introduce an upper limit on the structural volume to 
obtain useful results. For compliance minimization problems without design dependent loading a volume constraint 
is necessary to prevent the trivial solution where the design domain is completely filled with material. A volume 
constraint is not required here, as adding weight to the wing can increase the compliance because of the design 
dependent loads. Adding weight will increase the aerodynamic loads, through an increase angle of attack for trim, 
and also increase the body force loads. An upper limit on a may be required to prevent wing stall, which is highly 
undesirable in flight. The DLM cannot model stalled flow, and this effect is not included in this work. 

C. Shape Sensitivity Analysis 

The shape derivative for the problem defined in Eq. (29) is derived using results from Ref. 10 and accounting for 
the influence of the aerostructural coupling and the trim condition. The compliance can be written in a continuous 
form as: 


C(Q,u) = Nf (f g - u\dQ +f (f c - u)dT N - a(u ) f (f a ■ u)dT N + f {Q l jU l ■ Uj)dT N (30) 

J ^ u Tv u Tv u Tv x 

where T N is the part of the boundary where the aerodynamic loads are applied (such as the upper surface of the 
wing) and Q iyj is the continuous tensor form of the aerodynamic stiffness matrix. The continuous form of the 
aerostructural equilibrium equation is: 

f n ( Ee ( u )' £ ( v )) dQ = N fJfg + fr (/c -v)dr N -a[u) f (f a -v)dr N +f (Q u v r ujjdT N 

** ^ ^ 1 N 1 \ ** In (31) 

+ ( (v ' Ee(u)n + u ' Ee(v)n)dT D 

J r n 


where E is the material property tensor, e{u) is the strain tensor, v is a virtual displacement, n is a unit vector normal 
to the boundary and T D is the part of the boundary where displacement boundary conditions are applied. 

The adjoint method is used to obtain the shape derivative of Eq. (30). First, a Lagrangian function is formed by 
adding the aerostructural equilibrium equation, Eq. (31), to the objective, Eq. (30), using a Lagrange multiplier, A: 


A(Q,v,A) = C(Q,v)-fjEe(v)-e(A))dQ + NfjA-f g jd£2 + J {A- f c )dV N 

-a(v)f (A-f a )dT N + f (Q iJ A i -Vj)dr N +f (A-Ee(v)n + v-Ee(A)n)dT D 

J Tv J Tv x J E) 


(32) 


The adjoint equation is derived by differentiating Eq. (32) with respect to v in the direction 6 and defining ( u , p) 
as a stationary point, where 6 is an arbitrary vector. After collecting similar terms: 


dA 

dv 


(Q,u,p),e) = 0 = N f (f g -o)dQ+ f {f a (u)-6)dT N --^l f (f a ■ (u + p))6dT N 
/ JQ' 7 J T n du J T n 

+ f (Q i j(u i + p i )-dj)dT N - f [Es(e)-e(p))dQ + f [p- Ee(d)n + 6 ■ Ee(p)n)dT D 


(33) 


where the chevron brackets denote a vector inner product. To complete the derivative in Eq. (33) the change in angle 
of attack due to a change in deformed wing shape is required. From Eq. (21): 

da ( u)/du = b / L a (34) 
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Using FEA to solve the adjoint equation, Eq. (33) can be written in the following matrix form: 

K T p = N-f g +f a (u) + [Q T -(b®f a )/L a ](u + p), u = p = 0 on T D (35) 

The shape derivative can now be computed by differentiating the Lagrangian, Eq. (32) with respect to the shape 
in the direction 0 , at the stationary point ( u , p). During optimization, the part of the boundary where boundary 
conditions and aerodynamic loads are applied is fixed, thus 0 = 0 on these parts and we only compute the shape 
derivative on the remaining free part of the boundary, r 0 : 

^(Q,u,p)(6)= f e-n{N-f a -(u + p)-Ee(u)-e(p))dr 0 -^- [ (. f a -(u + p))dT N (36) 

oil dil •'Tv 

The second part of Eq. (36) is the change in the compliance with respect to a change in the angle of attack due to 
a change in the total weight of the structure. Thus, from Eq. (21) the derivative of a with respect to the shape, at the 
stationary point, is simply the shape derivative of the structural weight multiplied by the load factor over L a : 

da(u)/dQ = -N/L a f 6-n(pg)dT 0 (37) 

** r o 


where p is the density of the material and g is acceleration due to gravity. 

Setting On equal to the normal velocity, where a positive velocity is an inward movement: On = -V n , the final 
shape derivative on the free boundary becomes: 

OC r Nr 

— (Q ,u,p)= I V n (Ee(u) -e(p) - N ■ f a ■ (u + p) - A ■ pg)dT 0 , A = — I (/„•(« + p))dT N (38) 

ail •'T L a •'Tv 

where A is a factor, constant over the entire boundary, that relates the change in compliance due to a change in the 
angle of attack, resulting from a change in the wing box weight. 

To solve the aerostructural optimization problem stated in Eq. (29) using the level set method, the shape 
derivative must be negative to ensure a decrease in the objective. Thus, the velocity function is defined such that the 
shape derivative, Eq.(38) is negative everywhere on the free boundary. The velocity function is simply defined as: 

V n = N • f g - (u + p) + A- pg - Es(u)s(p) (39) 

This velocity function is used with the level set method described in Section II to progressively improve the 
structure with respect to the optimization objective. The sensitivity calculation was validated against sensitivities 
obtained using finite differences. 

D. Solution of the Primary and Adjoint Aerostructural Problems 

The aerostructural equilibrium equation, Eq.(29), can be restated as a linear equation and solved directly. 
However, in practice this is not an efficient approach, as the aerodynamic stiffness matrix is dense, whereas the 
structural stiffness matrix is sparse. Thus, methods usually employed to efficiently solve static linear elastic 
equations may not be efficient for solving the coupled equation, due to the dense matrix. Therefore, it can be more 
efficient to use an iterative approach to solve the aerostructural equilibrium equation. 

The iterative approach employed in this work starts with a load vector computed assuming no coupling between 
the aerodynamic loads and deformed wing shape. The deformed wing shape is obtained for the current load vector, 
which is then used to update the load vector by including the coupling terms. This process is repeated until the 
difference between the previous and updated load vectors is very small. The process can be summarized as: 

1 . Initialize load vector: / 0 = N-f g + f a ( 0) 

2. Solve: Ku k+1 =f k 

3. Update load vector: / k+1 = N-f g +f a (u k+l ) 

4. Check for convergence: ll/ k+1 -/ k II < small 
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The iterative approach is efficient, although it can be unstable and fail to converge, especially for more flexible 
wings. To overcome this, a not convergent solution is detected by tracking the relative difference in subsequent load 
vectors. If the solution is not converging, an under-relaxation approach is used to obtain the next load vector: / k+1 = 
N'f g + 0.25 f a (u k ) + 0.75 f a {u k ' 1 ). The relaxation factor of 0.25 was chosen to give a good compromise between 
efficiency and stability for the examples studied in this work. 

The iterative approach is also used to solve the adjoint equation, Eq.(35). In this case the displacement vector, u , 
is known from the solution of the primary equation. The initial adjoint load vector is computed using: p = 0. The 
adjoint solution is then obtained using the same iterative approach that is used to find the primary solution, including 
the under-relaxation method for stabilization. 

IV. Optimization of a CRM Wing Box 

The 3D unstructured grid level set method, detailed in Section II, is used solve the aero structural optimization 
problem, detailed in Section III. The subject of the optimization is the internal configuration of the CRM wing box. 28 



A. The Common Research Model Wing 

The CRM is a transonic transport configuration with an aspect ratio of 9, a taper ratio of 0.275, a sweep angle of 
35°, and a cruise Mach number of 0.85. For all flight loads a dynamic pressure of 5,897 Pa is used. The wing box is 
located between 12% and 71% of the local airfoil chord. The fixed weight of the aircraft is 500 kN, thus, in steady 
level flight, each wing must produce enough lift to support 250 kN plus the weight of one wing. The wing box 
structure is made from aluminum, with a Young’s modulus of 68.95 GPa, Poisson’s ratio of 0.3 and density of 2800 
kg/m 3 . 

The CRM wing box mesh is shown in Fig. 4 and has 150,300 hexahedral elements and 166,160 nodes. There are 
15 elements though the depth, 30 along the chord and 334 along the span. This mesh is used for the FEA and to 
discretize the implicit function. All degrees of freedom for nodes at the wing root are fixed. The DEM mesh is 
composed of 20 x 100 boxes in the chord- wise and span- wise directions, respectively, Fig. 5. The FPS method is 
used to apply the aerodynamic loads to a subset of 1705 nodes on the upper surface of the structural mesh, Fig. 6. 
When the FPS method is used to transform displacements into downwash, information is used from both the upper 
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and lower surfaces of the wing box mesh. The complete aeroelastic model, described in Section III, was validated 
against results obtained from NASTRAN. FEA and DLM results were validated independently, before validating the 
whole model, including the FPS method. 

The hexahedra finite elements used to discretize the wing box cannot model thin, shell-like structures typical of 
the conventional semi-monocoque configuration. Therefore, the structures obtained using the mesh in Fig. 4 are 
likely to be much heavier than real wing boxes. This will result in a high angle of attack and large body force loads, 
which are undesirable and make comparison to conventional configurations difficult. To obtain reasonable 
topologies with the mesh in Fig. 4, the modulus and density of the material used during the entire optimization are 
10% of the real values. 




Figure 5. DLM mesh (CRM wing box outlined in red). 


Figure 6. Down selected nodes used in the FPS method. 


The 3D level set optimization method is free to choose any structure within the wing box domain and is not 
restricted by any predefined configuration. However, during the optimization the outer layer of elements is fixed to 
remain part of the structure, therefore the upper and lower layers of elements act as the upper and lower skins. This 
ensures that holes do not develop in the upper and lower skins and that the skins remain connected. Finally the 
convergence criterion used for all examples is y = 8x1 O' 4 , Eq. (15). This value was chosen based on the observed 
convergent behavior of the examples studied in this work. 

B. Single Load Case 

This section considers several single load case problems for a range of flight maneuvers, characterized by the 
load factor, N , and a ground load case, which is a 2.5g taxi-bump. The ground load case only considers static body 
force loads and does not consider any aerodynamic loading or the trim condition. The initial configuration in each 
case is shown in Fig. 7, where material fills the wing box with regularly distributed spherical voids. The initial voids 
are required, as the current method cannot automatically create new holes during optimization. 

7. Steady level flight 

The first single load case problem is 
steady level flight with a load factor: N = 

1. The solution is shown in Fig. 8. The 
topology of the solution is relatively 
simple, as material is added to the 
outboard region and removed from the 
inboard region, although there are some 
subtle variations of the skin thickness in 
the inboard area. The extra material 
towards the tip helps reduce compliance 
by inertia relief, as the increase in the 
body force loading outboard acts to 
oppose the wing bending caused by the 
aerodynamic lift force. Although optimal 
for flight, the concentration of material at 

13 

American Institute of Aeronautics and Astronautics 







Figure 7. CRM wing box initial design. 


the tip is undesirable when considering ground load cases, such as the taxi-bump case, as the moment due to the tip 
mass can exceed that caused by the flight loads. 40 The material removed from the inboard region also helps reduce 
compliance, as the wing box weight is reduced, which results in a reduced angle of attack and reduced aerodynamic 
loading, Fig. 9. The two primary factors driving the topology in this case appear to be an increase of body force 
loading outboard and an overall weight reduction to reduce the root angle of attack and associated loading. The 
results show no significant aeroelastic tailoring during the optimization, as the ratio of tip deflection to tip twist 
changes by only 2.8% from the initial to the final structure. 



The unstructured 3D level set method performed well and the solution converged in 48 iterations, Fig. 10. 



Figure 9. Steady level flight, weight and angle of attack. Figure 10. Steady level flight, compliance. 


2. Flight maneuvers 

Two flight maneuvers are considered, a pull-up maneuver, with N = 2.5, and a push-over maneuver, with N = -1. 
These were chosen as they are typically used to size an aircraft at the boundaries of the flight envelope. 41 All other 
parameters are the same as for the steady level flight example. The solutions for the pull-up and push-over 
maneuvers are shown in Fig. 11 and 12, respectively. The topologies for the maneuver cases are very similar to the 
solution for the steady level flight case, Fig. 8, with only small differences in the thickness of the skins. This is not 
surprising, as the two main factors that drive the topology are the body force loading and aerodynamic loading from 
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the root angle of attack, both of which scale linearly with the load factor. The results are not exactly the same, as the 
constant part of the aerodynamic load vector does not vary with the load factor. It can be seen from the convergence 
histories of wing box weight, angle of attack and compliance, that all flight load cases follow a similar path during 
optimization, Figs 9, 10, 13-16, with only a noticeable difference in the magnitudes (and sign for the push-over 
maneuver angle of attack). 



Figure 11. Pull-up maneuver solution. 


Figure 12. Push-over maneuver solution. 



0 10 20 | lerarton 30 40 50 

Figure 13. Pull-up, weight and angle of attack. 


10 20 Ik-mtiun 30 ■W 

Figure 14. Pull-up, compliance. 


50 


3. Taxi-bump 

The load factor for the taxi-bump case is N = -2.5. No aerodynamic loads are applied and the trim condition is 
not enforced. The taxi-bump case is self-adjoint, as the adjoint equation, Eq. (35), reduces to the static equilibrium 
equation, thus u - p. The solution for the taxi-bump case is shown in Fig. 17 and is very different from the flight 
load case solutions. Material has been removed from the outboard region to reduce the inertia effect on the bending 
and some material has been added to the inboard region, to reinforce the wing box against bending. Some initial 
voids in the inboard region persist, which suggests that the removal of these voids has little impact on the 
compliance. Another difference from the flight load solutions is the lack of local variations in the wing skin 
thickness, which suggests that this feature is a result of the aeroelastic coupling for the flight load examples. Note 
that the solution shown in Fig. 17 has not reached the specified convergence criterion and was terminated after 400 
iterations, although the total compliance is not changing significantly at this point, Fig. 18. 
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C. Multiple Load Cases 

The results in the previous section show that optimal designs for flight load cases have a concentration of 
material in the outboard region for inertia relief, Figs. 8 , 11 and 12. This feature is undesirable for ground load cases 
and the optimal design for the taxi-bump case has minimal material outboard, Fig. 17. Therefore, the optimal 
designs for flight load and ground load cases are conflicting. The trade-off between these two load cases is explored 
by solving a multiple load case problem where one load case is steady level flight and the second load case is the 
taxi-bump: 


Minimize : rjCf + (l - r])C g (40) 

where C/ is the compliance from the steady level flight case, C g is the compliance from the ground load taxi-bump 
case and rj is a weighting factor for the two cases, where rj E [0,1]. The shape derivative for the multiple load case 
problem is simply the weighted sum of the derivatives from each case. 

The multiple load case problem, Eq. (40), was solved for five values of 77 : 0.0, 0.25, 0.5, 0.75, 1.0. The solution 
for 77 = 1.0 is the same as the steady level flight solution, Fig. 8 , and the solution for rj = 0.0 is the same as the taxi- 
bump solution, Fig. 17. Solutions for the other 77 values are shown in Fig. 19. Table 1 shows a summary of wing box 
weight values and compliance values, for each solution and both load cases, as expected the compliance for the 
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flight load case increases with decreasing rj and the opposite occurs for the ground load compliance. However, the 
trend for wing box weight is not so clear. 

Table 1. Multiple load cases, compliance and wing box weight values. 


V 

C f , kNm 

C s , kNm 

W„,kN 

1.00 

23.0 

39.8 

93.8 

0.75 

27.6 

9.0 

77.5 

0.50 

44.0 

5.2 

161.8 

0.25 

46.4 

3.5 

153.6 

0.00 

50.4 

3.2 

151.9 



Figure 19. Multiple load case solutions: a) r] = 0.75, b) rj = 0.5, c) rj = 0.25. 


The solutions for rj = 0.5 and 0.25 are similar to the taxi-bump solution, as material is removed from the 
outboard region. More material is removed as rj is decreased, as the inertial effects of the taxi-bump load case 
influence the optimization more. There are also some subtle skin thickness variations that are not seen for the pure 
taxi-bump solution, which may be a result of the aeroelastic coupling in the optimization. The compliance values for 
the rj = 0.5 and 0.25 solutions are within 10% of those for the taxi-bump case, where the flight load compliance is an 
order of magnitude greater than the ground load compliance, Table 1. The compliance convergence history for the rj 
= 0.5 problem shows that the flight load compliance actually increases during optimization, Fig. 20. The results 
suggests that the rj = 0.5 and 0.25 weighting values do not provide a good compromise between the flight and 
ground load cases for compliance minimization, as the solutions are not significantly different from the taxi-bump 
case and compliance values for the two load cases are an order of magnitude different. 


so 




0 25 50 75 JQ0 125 350 

I tiTaliim 


Figure 21. rj = 0.25, compliance. 
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The solution for 77 = 0.75 is significantly different to all previous solutions, as material is removed along the 
entire span of the wing, except for the tip leading edge region, where material is retained. Also, the compliance for 
both load cases is reduced during the optimization, Fig. 21, in contrast to the 77 = 0.5 result, Fig. 20. A detailed view 
of the solution for 77 = 0.75 that uses cross sections to show the material distribution of the interior structure is 
shown in Fig. 22. The skin thickness near the trailing edge reduces in the outboard region, whereas a significant 
amount of material is retained in the tip leading edge region. The extra material at the tip provides inertial relief for 
the flight load case, although less material is added compared with the solution for pure flight loads (Fig. 8 ) as tip 
mass is penalized by the taxi-bump load case. The increase of material towards the leading edge suggests that the 
optimization is finding advantage from an aeroelastic effect, as the elastic axis is moved forwards. This increases the 
natural wash-out of the back swept wing, leading to a reduced angle of attack outboard and hence the load 
distribution is moved inboard. This is highlighted in Fig. 23, which shows the DLM box pressure coefficients for the 
steady level flight load case for the 77 = 1.0 and 77 = 0.75 solutions. The pressure at the tip leading edge region is 
lower for the 77 = 0.75 solution, compared with the 77 = 1.0 solution. Moving the load distribution inboard is 
desirable for compliance optimization and wing structure design in general, as the overall bending moment on the 
wing is reduced. 

The weight of the 77 = 0.75 design is the lowest of all the optimal designs obtained in this work and the 
compliance values for the two load cases are much closer, with the flight load compliance approximately three times 
the ground load compliance, Table 1. Also, the fight load compliance is only 20% greater, but the taxi-bump 
compliance is reduced by 77% compared with the pure flight load optimum. The results suggest that the design 
obtained using 77 = 0.75 provides a good compromise between the flight and ground load cases for compliance 
minimization. 



Figure 22. Cross sections through the wing box (not ribs), r\ - 0.75 solution. 


D. Discussion 

The wing box structures obtained in this work using a 3D level set topology optimization method are very 
different from conventional structures. There are no spar or rib type members and the topologies mainly consist of 
large volumes that are either void or filled with material. There are several possible explanations for why solutions 
of this type are obtained with the method presented in this work. First, it appears that the solutions are often driven 
by the inertia effects and, for the flight load cases, the aerodynamic loads associated with the root angle of attack for 
trim. The wing is cantilevered from its root and, therefore, the influence of the body force loading is likely to drive 
the optimizer to simply add or remove material in the outboard region, where inertia effects are more significant. 
The angle of attack for trim is mainly influenced by the weight of the wing box, which is a global quantity and 
therefore unlikely to produce local features during optimization. However, the analysis and shape derivative do 
account for the aeroelastic feedback on the structural deformation, angle of attack and, therefore, the total 
compliance. The influence of the aeroelastic feedback on the optimized solutions for the examples studied in this 
work often appears less significant than the influence of the body force loading and total wing box weight. However, 
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evidence of the aeroelastic feedback may be seen in the variations in the skin thickness for the flight load case 
examples and especially for the multiple load case problem with rj = 0.75. Further work is required to understand the 
influence of each factor on the optimized designs. It would also be interesting to include other loads, such as weight 
from fuel, engines and leading and trailing edge devices and additional aerodynamic loads from flaps and ailerons. 

The role of ribs and spars in conventional wing box structures is not simply to provide stiffness and there are 
many constraints, such as stresses, buckling, aeroelastic stability and maintaining the aerodynamic shape of the 
wing. The problem studied in this work does not take into account such constraints that could drive the structure to 
form more spar/rib-like internal members. 

A further consideration is the part of the wing box structure that is not part of the design space. The outer layer 
of elements was forced to remain part of the structure during the entire optimization, primarily to provide the 
aerodynamic surface. The mesh was reasonably coarse; therefore, the contribution of the fixed elements to the 
overall stiffness is significant and could be influencing the optimal design of the wing box interior. This factor could 
be investigated by using a finer mesh, so that the influence of the fixed outer layer of elements is reduced. Another 
approach could be to include the outer layer, or at least the wing skins, as part of the design space, although the skins 
should remain continuous, without holes, to support the aerodynamic loads. 

The solutions obtained may also be influenced by the limitations of the level set based optimization method. The 
method employed in this work does not have a mechanism for introducing new holes. Therefore, the solutions are 
dependent on the initial design. 11 ' 13 This factor could be investigated by starting from different initial designs or by 
introducing hole creation. 13 



Figure 23. DLM box pressure coefficients: a) steady level flight solution, b) rj = 0.75 solution. 


V. Conclusions 

This paper uses level set based topology optimization to design the wing box structure of the Common Research 
Model. The optimizer is free to choose any structure within the 3D wing box space and the level set method is 
employed on a 3D unstructured grid. A fast marching method is employed and developed to efficiently compute the 
velocity extension and implicit function reinitialization on the unstructured grid. 

The compliance of the wing box is minimized under body force and aerodynamic loading, which is coupled to 
the deformed wing shape. The root angle of attack is allowed to change during the optimization to meet the trim 
condition that weight equals lift. The adjoint method is used to obtain the shape derivative of the coupled 
aero structural problem, which is then used to define the velocity function used in the level set method. 

The compliance minimization problem is solved for single load cases, including steady level flight, maneuvers 
and a ground load tax-bump case. The solutions for the flight load cases are very similar, with material added to the 
outboard region for inertia relief. In contrast, the solution for the taxi-bump case removed material from the 
outboard region to reduce the inertial loading effect. To resolve this conflict in optimal deigns for flight and ground 
load cases, a multiple load case problem is formulated, where one case is steady level flight and the other is the taxi- 
bump. A load case weight factor of 0.75 (on the flight load case) provided a reasonable compromise solution for 
good flight and ground load performance. 
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The method is able to efficiently and smoothly obtain solutions for the compliance minimization of the CRM 
wing box. The solutions are remarkably different from conventional rib- spar type designs, as the topologies mainly 
consist of large volumes that are either void or filled with material. Further refinements are suggested to investigate 
this, such as, the problem definition to include constraints, the modeling of the fixed part of the part of the structure 
and the lack of a hole creation mechanism in the level set method. 
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